/*
 * MyLittleDam.cpp
 *
 *  Created on: 26/ott/2012
 *      Author: lordmzn
 */

#include "MyLittleDamSimulator.h"
#include "moeaframework.h"
#include <iostream>
#include <cmath>
#include <cstddef>
#include <vector>
#include <algorithm>
#include <random>

#ifdef RANDOM_NUM_100
static double numCasuali[100] = { 41.554336, 34.782222, 39.449927, 38.729340,
		41.819383, 42.046036, 43.011905, 44.009569, 39.871777, 37.192379,
		41.542659, 37.738619, 45.874758, 44.962047, 40.648001, 45.256974,
		38.268144, 34.182198, 43.225501, 38.906461, 35.772933, 37.978430,
		38.761098, 40.285301, 38.764016, 42.671219, 35.176639, 34.070290,
		43.592754, 37.896099, 41.460852, 36.392245, 34.263731, 38.930609,
		44.163730, 41.583630, 42.078223, 37.919533, 37.339892, 40.201364,
		39.265638, 40.310803, 40.416767, 39.192028, 42.862306, 37.890663,
		42.163035, 40.898607, 44.614258, 41.845280, 34.182367, 37.022532,
		40.882168, 38.319578, 39.911324, 40.404525, 41.961012, 40.197465,
		37.846977, 40.423441, 38.131560, 37.620548, 43.241820, 47.091243,
		38.199274, 39.944762, 44.552388, 38.173059, 40.423040, 37.551460,
		36.292238, 38.214299, 45.078573, 33.927022, 45.500637, 36.889676,
		41.086988, 35.276466, 45.716899, 40.234793, 31.212507, 42.959565,
		42.471606, 40.854856, 38.580483, 38.122061, 41.119499, 38.098986,
		39.965605, 40.776878, 43.816923, 34.184480, 36.758442, 41.411816,
		41.773195, 42.517080, 40.814041, 41.242223, 43.424939, 36.902555 };
#endif
#ifdef RANDOM_NUM_200
static double numCasuali2[200] = {42.348853, 40.385337, 40.263645, 39.625107,
	39.190262, 39.009182, 42.593639, 36.570410, 33.811144, 35.902019,
	40.973346, 39.201676, 43.079371, 42.810063, 41.719289, 37.871780,
	40.673814, 39.189500, 36.390448, 39.642664, 44.721819, 40.811457,
	32.852510, 37.859464, 38.354802, 38.338110, 42.066871, 40.713699,
	33.440446, 41.010544, 40.149490, 41.033329, 40.827867, 41.485764,
	39.255500, 41.907416, 35.542261, 35.986264, 36.617384, 40.856662,
	40.154229, 46.992240, 33.883007, 39.380989, 36.113300, 40.569301,
	41.583800, 41.101421, 40.180565, 39.022753, 42.769191, 35.511317,
	40.668923, 42.337482, 37.332473, 39.902737, 31.687839, 45.737732,
	45.671909, 41.334343, 40.274383, 37.995208, 40.724667, 42.181605,
	42.237647, 40.822237, 37.584129, 43.661456, 45.724612, 39.374550,
	41.843862, 35.102335, 38.849194, 41.203342, 35.246264, 39.705061,
	38.178010, 37.144202, 36.863201, 42.907110, 41.130906, 38.497772,
	38.575600, 35.127518, 42.597008, 38.064684, 35.237799, 44.253072,
	36.786134, 44.807414, 31.783756, 43.422833, 37.023345, 41.646825,
	42.237966, 33.284303, 33.850668, 38.019926, 38.915011, 39.576424,
	37.969102, 37.586738, 36.867017, 42.509949, 39.616618, 38.263951,
	47.396929, 36.315048, 43.714778, 43.548247, 39.717323, 36.132280,
	43.082552, 37.221569, 42.059050, 38.607678, 35.860941, 45.672841,
	40.888675, 41.283308, 37.504841, 43.226186, 38.691598, 39.360953,
	41.196816, 40.724447, 37.824073, 40.405706, 39.864480, 42.100969,
	40.679967, 42.152113, 39.460325, 41.666326, 42.156830, 37.237522,
	41.602952, 36.450688, 38.185030, 37.316167, 40.912332, 39.181667,
	39.933748, 43.911276, 37.687809, 44.608864, 40.607173, 36.265507,
	39.624198, 40.546378, 34.650808, 40.147035, 38.197810, 38.286045,
	38.643790, 37.683719, 41.692905, 45.781326, 42.126263, 37.753258,
	40.379056, 42.295546, 39.330776, 39.905880, 37.420365, 38.992640,
	37.513334, 42.037389, 40.881602, 41.099287, 39.856077, 45.082347,
	41.519656, 37.043604, 36.019611, 42.400070, 36.268647, 39.256020,
	35.986122, 41.979485, 41.590744, 41.966787, 36.712334, 39.134395,
	41.926683, 38.666782, 37.411425, 36.961151, 41.293573, 40.802289,
	41.194899, 36.781833, 45.290942, 42.241662, 45.346207, 44.904726,
	41.117299, 38.574522, 42.862424, 38.426358};
#endif

MyLittleDamSimulator::MyLittleDamSimulator(double mean, double stddev,
		int length) :
#ifdef RANDOM_NUM_100
				inflow(numCasuali,
						numCasuali + sizeof(numCasuali) / sizeof(numCasuali[0])) {
#endif
#ifdef RANDOM_NUM_200
	inflow(numCasuali2,
			numCasuali2 + sizeof(numCasuali2) / sizeof(numCasuali2[0])) {
#endif
#ifdef RANDOM_GENERATOR_MONTECARLO_10
	inflow(length, 0.0) {
		extractAfflussi(mean, stddev);
#endif
}

MyLittleDamSimulator::~MyLittleDamSimulator() {

}

void MyLittleDamSimulator::extractAfflussi(double mean, double stddev) {
	std::normal_distribution<double> normal(mean, std::sqrt(stddev));
	std::default_random_engine generator(12);
	for (unsigned int i = 0; i < inflow.size(); i++) {
		inflow[i] = normal(generator);
	}
}

double MyLittleDamSimulator::calcolaDecisione(double p1, double storageValue) {
	return (double) p1 * storageValue;
}

double MyLittleDamSimulator::RBFs_policy(std::vector<double> input,
		std::vector<double> parameters, unsigned int N) {
	// setting
	unsigned int D = input.size(); // number of states
	const double u_min = 0.0;
	const double u_max = 100;
	// reshape of parameters
	double c[N][D]; // center = (N-RBFs, D-state variables)
	double b[N][D]; // radius = (N-RBFs, D-state variables)
	unsigned int count = 0;
	for (unsigned int i = 0; i < N; i++) {
		for (unsigned int j = 0; j < D; j++) {
			c[i][j] = parameters[count];
			b[i][j] = parameters[count + 1];
			count = count + 2;
		}
	}
	std::vector<double> w; // weights of the RBFs
	double ws = 0;
	for (unsigned int i = 0; i < N - 1; i++) {
		w.push_back(parameters[count]);
		ws = ws + parameters[count];
		count = count + 1;
	}
	w.push_back(1 - ws);
	// Normalization of the input
	double input_norm[D][2];
	input_norm[0][0] = 0;
	input_norm[0][1] = 200;
	//input_norm[1][0] = 0 ;
	//input_norm[1][1] = 200 ;
	std::vector<double> input1;
	for (unsigned int i = 0; i < D; i++) {
		input1.push_back(
				(input[i] - input_norm[i][0])
						/ (input_norm[i][1] - input_norm[i][0]));
	}
	// RBFs
	std::vector<double> phi;
	double bf = 0;
	for (unsigned int j = 0; j < N; j++) {
		for (unsigned int i = 0; i < D; i++) {
			bf = bf
					+ (input1[i] - c[j][i]) * (input1[i] - c[j][i])
							/ (b[j][i] * b[j][i]);
		}
		phi.push_back(exp(-bf));
	}
	// compute decision
	double out = 0;
	for (unsigned int i = 0; i < N; i++) {
		out = out + w[i] * phi[i];
	}
	return out * (u_max - u_min) + u_min;
}

void MyLittleDamSimulator::runSimulation(double vars[], double objs[],
		DamSimParam* p) {
#ifdef RANDOM_GENERATOR_MONTECARLO_10
	double obj0 = 0.0;
	double obj1 = 0.0;
	for (int i=0; i < MONTECARLO_NRUNS; i++ ) {
#endif
	std::vector<double> storage(p->horizon + 1, 0.0);
	std::vector<double> irr_stepcost(p->horizon + 1, 0.0);
	std::vector<double> eso_stepcost(p->horizon + 1, 0.0);
	std::vector<double> params(p->nvars, 0.0);
	storage[0] = p->s0;
	double ril_SUP = 0.0;
	double ril_INF = 0.0;
	double release;
	double decision = 0.0;
	for (int i = 0; i < p->nvars; i++) {
		params[i] = vars[i];
	}

	for (std::vector<double>::size_type t = 0; t < p->horizon; t++) {
		//std::cout << "Time step: t \t\t= " << t << ";" << std::endl;
		//std::cout << "State (storage): s(t) \t= " << storage[t] << ";" << std::endl;
		decision = RBFs_policy(std::vector<double>(1, storage[t]), params,
				p->RBF_N);
		//std::cout << "Control: u(t+1) \t= " << decisione << ";" << std::endl;
		ril_SUP = storage[t];
		ril_INF = std::max(storage[t] - 100, 0.0);
		release = std::min(ril_SUP, std::max(ril_INF, decision));
		//std::cout << "Release r(t+1) \t= " << release << "; R_INF = " << ril_INF << "; R_SUP = " << ril_SUP << ";" << std::endl;
		//std::cout << "Afflusso: a(t+1) \t= " << afflusso[t+1] << ";" << std::endl;
		storage[t + 1] = storage[t] + inflow[t + 1] - release;
		irr_stepcost[t + 1] = std::pow(p->water_request - release, 2);
		eso_stepcost[t + 1] = std::max((storage[t + 1] - p->liv_crit), 0.0);
	}

	objs[0] = 0.0;
	objs[1] = 0.0;
	for (int t = 0; t < (int) irr_stepcost.size(); t++) {
		objs[0] = objs[0] + irr_stepcost[t];
		objs[1] = objs[1] + eso_stepcost[t];
	}
	objs[0] = objs[0] / (double) irr_stepcost.size();
	objs[1] = objs[1] / (double) eso_stepcost.size();
#ifdef RANDOM_GENERATOR_MONTECARLO_10
	obj0 = obj0 + objs[0];
	obj1 = obj1 + objs[1];
	extractAfflussi(p->aff_mean, p->aff_stddev);
}
objs[0] = obj0 / MONTECARLO_NRUNS;
objs[1] = obj1 / MONTECARLO_NRUNS;
#endif
}

int main(int argc, char **argv) {
	DamSimParam* param = new DamSimParam();
	param->horizon = 100;
#ifdef RANDOM_NUM_200
	param->horizon = 200;
#endif
	param->s0 = 50.0;
	param->water_request = 50.0;
	param->liv_crit = 50.0;
	param->aff_mean = 40.0;
	param->aff_stddev = 10.0;
	param->RBF_N = 2;
	param->nvars = 5;
	param->nobjs = 2;
	double vars[param->nvars];
	double objs[param->nobjs];

	MyLittleDamSimulator* dam = new MyLittleDamSimulator(param->aff_mean,
			param->aff_stddev, param->horizon);

	MOEA_Init(param->nobjs, 0);

	while (MOEA_Next_solution() == MOEA_SUCCESS) {
		MOEA_Read_doubles(param->nvars, vars);
		dam->runSimulation(vars, objs, param);
		MOEA_Write(objs, NULL);
	}

	MOEA_Terminate();

	delete param;
	delete dam;

	return 0; // was: return EXIT_SUCCESS;
}

